/*******************************************************************************

        SSC0143 - Programacao Concorrente
        Prof. Julio Cesar Estrela
       
        T1 - Calculo de PI com 10 000 000 de casas decimais
       
        Grupo:
       
        Daniel de Oliveira Pereira              6427171
        Guilherme Takasawa Yagui                6426698
        Vinicius Marques Stocco                 6427233
       
        Março de 2013
       
********************************************************************************        
                 Algoritmo de Gauss-Lagendre - Versao Concorrente      
*******************************************************************************/

#include <stdio.h>
#include <math.h>
#include <time.h>
#include <gmp.h>
#include <pthread.h>
#include <semaphore.h>

// numero de casas decimais corretas de PI
#define PREC_DGT        10000000

// numero de bits minimo para armazenar PI com o numero de casas decimais corretas
#define PREC_BIT        (int) ceil((log(10) / log(2)) * PREC_DGT)

// numero de iteracoes necessario para calcular PI com o numero de casas decimais corretas desejado
#define ITR             (int) ceil(log(PREC_DGT) / log(4))

mpf_t an, an1, bn, bn1, tn, tn1, pn, pn1, pi; //calculo pi
mpf_t tmp_a, tmp_b, tmp_t, tmp_p, aux; // auxiliares

sem_t s_main, s_A, s_B, s_T, s_P;   //semaforos

pthread_t t_a, t_b, t_t, t_p; //threads

// calcula a(n+1)
// a(n+1) = (a(n) + b(n))/2
void *A(void *args)
{      
        while (1)
        {
                // auto-bloqueio
                sem_wait(&s_A);

                mpf_add(tmp_a, an, bn);		
		mpf_div_ui(tmp_a, tmp_a, 2);	
       
                // libera a thread principal
                sem_post(&s_main);
        }
}

// calcula b(n+1)
// b(n+1) = sqrt(a(n)*b(n))
void *B(void *args)
{
        while (1)
        {
                // auto-bloqueio
                sem_wait(&s_B);
       
                mpf_mul(tmp_b, an, bn);
		mpf_sqrt(tmp_b, tmp_b);
       
                // libera a thread principal
                sem_post(&s_main);
        }
	
}


// calcula t(n+1)
// t(n+1) = t(n) - p(n)*(a(n) - a(n+1))^2
void *T(void *args)
{
	while (1)
        {
                // auto-bloqueio
                sem_wait(&s_T);

                mpf_sub(tmp_t, an, an1);
		mpf_pow_ui(tmp_t, tmp_t, 2);
		mpf_mul(tmp_t, pn, tmp_t);
		mpf_sub(tmp_t, tn, tmp_t);
       
                // libera a thread principal
                sem_post(&s_main);
        }
	
}

// calcula p(n+1)
// p(n+1) = 2*p(n)
void *P(void *args)
{
	while (1)
        {
                // auto-bloqueio
                sem_wait(&s_P);
       
                mpf_mul_ui(tmp_p, pn, 2);
       
                // libera a thread principal
                sem_post(&s_main);
        }
	
}


int main(char argc, char *argv[])
{
        clock_t t;              // tempo
	int i, f;
       
        // define a precisão desejada para cada variavel
        mpf_set_default_prec(PREC_BIT);
       
        // marca  o clock inicial da execucao do algoritmo
        t = clock();
       
        // aloca memoria para as variaveis utilizadas para calcular PI
        mpf_init(an);
        mpf_init(an1);
        mpf_init(bn);
        mpf_init(bn1);
        mpf_init(tn);
        mpf_init(tn1);
        mpf_init(pn);
        mpf_init(pn1);
	mpf_init(aux);
        mpf_init(pi);
        mpf_init(tmp_a);
        mpf_init(tmp_b);
        mpf_init(tmp_t);
        mpf_init(tmp_p);

        // inicializa a(0) = 1
        mpf_set_d(an, 1.0);

        // inicializa b(0) = 1/sqrt(2)
	mpf_set_d(bn, 1.0);
	mpf_set_d(aux, 2.0);
	mpf_sqrt(aux, aux);
	mpf_div(bn, bn, aux);

	// inicializa t(0) = 1/4
	mpf_set_d(tn, 1.0);
	mpf_div_ui(tn, tn, 4);

        // inicializa p(0) = 1
	mpf_set_d(pn, 1.0);

	// inicializa os semaforos
        // libera main, A, B, P e bloqueia T
        sem_init(&s_main, 0, 0);
        sem_init(&s_A, 0, 1);
        sem_init(&s_B, 0, 1);
        sem_init(&s_P, 0, 1);
        sem_init(&s_T, 0, 0);

        // cria as threads utilizadas para calcular a(n+1), b(n+1), t(n+1) e p(n+1)
        pthread_create(&t_a, NULL, &A, NULL);
        pthread_create(&t_b, NULL, &B, NULL);
        pthread_create(&t_t, NULL, &T, NULL);
        pthread_create(&t_p, NULL, &P, NULL);

        for (i = 0; i < ITR; i++)
        {
		// bloqueia a thread main
                sem_wait(&s_main);
		sem_wait(&s_main);
		
                // calcula pi = ((a(n) + b(n))^2)/(4*t(n))
                mpf_add(pi, an, bn);
		mpf_pow_ui(pi, pi, 2);
		mpf_div_ui(pi, pi, 4);
		mpf_div(pi, pi, tn);

                gmp_printf("Iteracao %d, pi = %.*Ff \n", i, PREC_DGT, pi);

		// calcula a(n+1)
                mpf_set(an1, tmp_a);
		// calcula b(n+1)
                mpf_set(bn1, tmp_b);
		// calcula p(n+1)
                mpf_set(pn1, tmp_p);

		// libera a thread T
                sem_post(&s_T);

		// bloqueia a thread main
                sem_wait(&s_main);
		sem_wait(&s_main);

                // calcula t(n+1)
                mpf_set(tn1, tmp_t);

		//a(n) = a(n+1)
		//b(n) = b(n+1)
		//t(n) = t(n+1)
		//p(n) = p(n+1)
                mpf_set(an, an1);
                mpf_set(bn, bn1);
                mpf_set(tn, tn1);
                mpf_set(pn, pn1);

                // libera as threads
                sem_post(&s_A);
                sem_post(&s_B); 
                sem_post(&s_P); 
        }

        printf("\nNumero de iteracoes : %d \n", ITR);
       
        // cancela a execucao das threads
        pthread_cancel(t_a);
        pthread_cancel(t_b);
        pthread_cancel(t_t);
        pthread_cancel(t_p);
       
        // libera recursos alocados pelas threads
        pthread_join(t_a, NULL);
        pthread_join(t_b, NULL);
        pthread_join(t_t, NULL);
        pthread_join(t_p, NULL);
       
        // destroi os semaforos
        sem_destroy(&s_main);
        sem_destroy(&s_A);
        sem_destroy(&s_B);
        sem_destroy(&s_T);
        sem_destroy(&s_P);
       
        // libera a memoria alocada para o calculo de PI
        mpf_clear(an);
        mpf_clear(bn);
        mpf_clear(tn);
        mpf_clear(pn);
        mpf_clear(pi);
        mpf_clear(tmp_a);
        mpf_clear(tmp_b);
        mpf_clear(tmp_t);
        mpf_clear(tmp_p);
       
        // calcula o numero clocks necessario para executar o algoritmo
        t = clock() - t;
       
        printf("\nTempo de execucao : %.3f s \n\n", (double) t / CLOCKS_PER_SEC);

        return 0;
}

